Outline

• Anouncements
• RNA-seq Analysis Workflow
  • Read Alignment
  • Differential Expression
• RNA-seq Tools
• Odyssey Fest

Anouncements

• Please do not share your code with fellow students via github, please make sure your repo is private.
• Please double check you have access to Odyssey and /n/stat115/. You will need it for homework 3, 4, and 5 (due 3/4). Please shoot me an email if you are not blessed with access to that.
• Start early early early on homework 3!

RNA-seq analysis workflow

Wang, Z., Gerstein, M., & Snyder, M. (2009). RNA-Seq: a revolutionary tool for transcriptomics. Nature reviews genetics, 10(1), 57.

Fragment alignment

• Star: ultrafast universal RNA-seq aligner
  • Suffix array (SA) approach to leverage speed of binary searching
  • Requires an index (the SA) to be created prior to aligning
• Salmon: fast and bias-aware quantification of transcript expression using dual-phase inference
  • Pseudo-alignment concept to rapidly align reads
  • Reads are classified according to their fragment equivalence class
  • (Optionally) removes sequencing biases and computes coverage via likelihood maximization
  • Similar in concept and result to Near-optimal probabilistic RNA-seq quantification

Star alignment

• Needs to be run on Odyssey HPC
• Takes several hours to complete (given an index)
• How do we do that again?

#!/bin/bash
#SBATCH -N 1                   # Number of nodes
#SBATCH -n 10                  # Number of cores
#SBATCH -t 180                 # Runtime in minutes (0~10080)
#SBATCH -p general             # Partition
#SBATCH --mem=50000            # Total memory (varies across nodes)
#SBATCH -o star_%j.out         # Standard out goes to this file
#SBATCH -e star_%j.err         # Standard err goes to this file
#SBATCH --mail-type=END        # Email
#SBATCH --mail-user=YOUR_EMAIL

module load gcc/4.8.2-fasrc01 STAR/2.5.0c-fasrc02

STAR --genomeDir $GENOME \
     --readFilesIn $FASTQ1 $FASTQ2 \
     --outFileNamePrefix $OUTDIR/ \
     --outSAMprimaryFlag AllBestScore \
     --outSAMtype BAM SortedByCoordinate \
     --runThreadN 10 \
     --alignEndsType EndToEnd 

sbatch STARalignment.sh

Salmon tutorial

Salmon tutorial 1

How to train your Salmon

Salmon Tutorial

Salmon tutorial

Salmon tutorial 2

Salmon tutorial

Salmon tutorial 3

Salmon tutorial

Salmon tutorial 4

Salmon tutorial

Salmon tutorial 5

Salmon tutorial

Salmon tutorial 6

Salmon tutorial

Salmon tutorial 7

Salmon tutorial

Salmon tutorial 8

Salmon tutorial

Salmon tutorial 9

Salmon tutorial

Salmon tutorial 10

Salmon tutorial

Salmon tutorial 11

Salmon tutorial

Enjoy your Salmon

Salmon alignment

• Much quicker but will still be run on Odyssey
• First step is to create an index, then align
• For HW3, I have prepared the transcriptome file as /n/stat115/HW3_2019/transcriptome/Homo_sapiens.GRCh38.cdna.all.fa on the Odyssey

#!/bin/bash
#SBATCH -N 1                   # Number of nodes
#SBATCH -n 10                  # Number of cores
#SBATCH -t 180                 # Runtime in minutes (0~10080)
#SBATCH -p general             # Partition
#SBATCH --mem=50000            # Total memory (varies across nodes)
#SBATCH -o star_%j.out         # Standard out goes to this file
#SBATCH -e star_%j.err         # Standard err goes to this file
#SBATCH --mail-type=END        # Email
#SBATCH --mail-user=YOUR_EMAIL

module load salmon

salmon index -t $TRANSCRIPTOME -i $INDEX

sbatch createSalmonIndex.sh

#!/bin/bash
#SBATCH -N 1                   # Number of nodes
#SBATCH -n 10                  # Number of cores
#SBATCH -t 180                 # Runtime in minutes (0~10080)
#SBATCH -p general             # Partition
#SBATCH --mem=50000            # Total memory (varies across nodes)
#SBATCH -o star_%j.out         # Standard out goes to this file
#SBATCH -e star_%j.err         # Standard err goes to this file
#SBATCH --mail-type=END        # Email
#SBATCH --mail-user=YOUR_EMAIL

module load salmon

salmon quant -i $INDEX \
             -l A \
             -1 $FASTQ/ENCFF500PDO_sub.fastq\
             -2 $FASTQ/ENCFF708KQE_sub.fastq \
             -o $OUT \
             --numBootstraps 100 \
             -p 10 \
             --gcBias

sbatch Salmonalignment.sh

Running your own alignment

• Use the commands given earlier to run your own alignment in Homework 3
• Make full use of all possible partitions on the Odyssey
• Do manage your files. Construct a file to your local folder. Batch submit your job to the cluster, outputting files to your local folder
• Upon completion, you may receive an email
• View the slurm- file and the log file (located in the output directory)
• Abundance estimates are given in the .sf file. Copy these to your local directory (e.g. scp or fileZilla) for downstream analysis (DESeq part)

Differential expression

• Given the abundance estimates, we now wish to determine which genes/transcripts are differentially expressed between condition.
• Just like in microarrays, we need to create a design matrix specifying what we wish to estimate
• We will be making use of the DESeq2 packages. Install this via bioconductor.

# Install required packages
source("https://bioconductor.org/biocLite.R")
biocLite("BiocUpgrade")
biocLite("DESeq2")
biocLite("tximport")
biocLite("EnsDb.Hsapiens.v86")
biocLite("EnsDb.Mmusculus.v79")
install.packages("rjson")

Differential expression

• Load the package.
• Create a data frame with columns as specified below containing the necessary information for the design matrix (sampleName, fileName, condition).

library(DESeq2)
files <- grep("sf",list.files("Data"),value=TRUE)
condition <- c("4oh", "4oh", "4oh", "ctrl", "ctrl", "ctrl")
names <- c("4oh1", "4oh2", "4oh3", "ctrl1", "ctrl2", "ctrl3")
sampleTable <- data.frame(sampleName = files, fileName = files, condition = condition)

Differential expression

• Utilizing the mouse transcriptome as below, extract the transcript id columns and the gene id columns.
• Might need to change the Ensembl based annotation package for mouse into the package for human, say change EnsDb.Mmusculus.v79 to EnsDb.Hsapiens.v86.

library(EnsDb.Mmusculus.v79)
txdf <- transcripts(EnsDb.Mmusculus.v79, return.type="DataFrame")
tx2gene <- as.data.frame(txdf[,c("tx_id", "gene_id")])

Differential expression

• Now we’re going to import the data and format the data for analysis.
• Remove rows with a total row count less than or equal to 1.
• This step might take a few minutes.
• Run the analysis!

library(tximport)
txi <- tximport(file.path("Data",files), type="salmon", ignoreTxVersion = TRUE, tx2gene = tx2gene)
dds <- DESeqDataSetFromTximport(txi,colData=sampleTable,design=~condition)
dds <- dds[rowSums(counts(dds)) > 1, ]
dds <- DESeq(dds)

Differential expression

• Extract the significantly DE genes at the .05 significance level.
• Remove incomplete observations (those with missing column data).
• Find those DE genes significantly up/down regulated with an adjusted p-value below .05. How many are there?

res <- results(dds, alpha = 0.05)
res <- res[complete.cases(res),]
res <- res[order(res$padj),]
upR <- res[(res$padj < 0.05) & (res$log2FoldChange > 0),]
downR <- res[(res$padj < 0.05) & (res$log2FoldChange < 0),]
nrow(upR)

## [1] 382

nrow(downR) 

## [1] 518

Visualizing results

• Which plot should we use to visualize the differential expression? Use it to visualize the results.
• Any other thoughts on how to further visualize the results?

plotMA(res)

absOrdered <- rbind(upR,downR)
absOrdered <- absOrdered[order(abs(absOrdered$log2FoldChange),decreasing = TRUE),]
mostvariable <- log2(txi$abundance[row.names(absOrdered),]+.0001)

library(gplots)
heatmap.2(mostvariable[1:100,],trace="none",col=greenred(10))

Odyssey Fest

• A summary for our experiences

Odyssey Fest

• How to make full use of all partitions

#!/bin/bash
#SBATCH -N 1                   # Number of nodes
#SBATCH -n 10                  # Number of cores
#SBATCH -t 240                 # Runtime in minutes (0~10080)
#SBATCH -p general,serial_requeue,shared  # Partition
#SBATCH --mem=50000            # Total memory (varies across nodes)
#SBATCH -o salmon_%j.out         # Standard out goes to this file
#SBATCH -e salmon_%j.err         # Standard err goes to this file
#SBATCH --mail-type=END        # Email
#SBATCH --mail-user=ytang@hsph.harvard.edu

module load salmon

salmon index -t $TRANSCRIPTOME -i $INDEX

Basic Operations Reminder (A Small Cheat Sheet)

• Log in the cluster: ssh command

Basic steps to access the cluster

• Running jobs interactively: srun command

• View files under current location: ls command
• Change path: cd command
• Create path: mkdir command

• Rename path: mv command

Useful commands for data management

Cluster Computing Reminder (A Small Cheat Sheet)

• Upload scripts to the Odyssey

Inquire path on the Odyssey

Upload scripts to the Odyssey

• View scripts on the Odyssey

View script list

View specific script

View specific script

• Submit jobs to the cluster: sbatch command
• Manage jobs on the cluster: sbatch command
• Kill jobs on the cluster: scancel command

Manage current jobs

• Locate the output and download: scp command

Download the output

Good Luck!

• Office hour untill 7:00 pm